1 load packages and define colors

pacman::p_load(rmarkdown, tidyverse, phyloseq, magrittr, janitor, microbiome, knitr, lubridate, naniar, readxl, ggplot2, ggpubr,rstatix, metadeconfoundR, microViz)

id_palette <- c("11"="#4E79A7FF", "13"="#A0CBE8FF", "15"="#F28E2BFF", "16"="#FFBE7DFF", 
                "17"="#59A14FFF", "24"="#8CD17DFF", "25"="#B6992DFF", "26"="#F1CE63FF", 
                "27"="#499894FF", "29"="#86BCB6FF", "32"="#E15759FF", "33"="#FF9D9AFF", 
                "35"="#79706EFF", "5"="#BAB0ACFF", "6"="#D37295FF", "9"="#FABFD2FF", 
                "21"="#B07AA1FF", "31"="#D4A6C8FF", "8"="black")

2 read data

df <- read_excel("/Users/rebecca/Documents/Forschung/IMMProveCF/2024_09_03_Mohr test Nacl calculations (Nacho).xlsx", sheet = "Results")
## New names:
## • `` -> `...6`
## • `` -> `...10`
## • `` -> `...14`
## • `` -> `...15`
## • `` -> `...16`
## • `` -> `...17`
## • `` -> `...18`
## • `` -> `...19`
## • `` -> `...20`
## • `` -> `...21`
## • `` -> `...22`
## • `` -> `...23`
## • `` -> `...24`
## • `` -> `...25`
## • `` -> `...26`
## • `` -> `...27`
ps_sputum_full <- readRDS("/Users/rebecca/Documents/Forschung/IMMProveCF/R_analysis/rds_files/ps_IMProveCF_patients_Run1-23_18102023.rds")

df <- clean_names(df)

df <- df %>% 
  filter(!is.na(eppi_n))

df <- df[2:54, ]

df <- df %>% 
  select(-c(na_cl_conc, x14)) %>% # those are meaningless
  rename(na_cl_conc_mmolL_TIT=x15) %>% 
  rename(comments=x16) %>% 
  rename(na_cl_conc_mmolL_ISE=x19) %>% 
  rename(na_cl_conc_mmolL_COND=x23) %>% 
  rename(confidence_titration=x26) %>% 
  rename(confidence_anlytics=x27)

df <- df %>%
  mutate(across(8:21, as.numeric)) %>% 
  mutate(visit = as_factor(gsub("V", "", sample))) %>% 
  mutate(visit = factor(visit, levels = sort(as.numeric(levels(factor(visit))))))
## Warning: There were 3 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `across(8:21, as.numeric)`.
## Caused by warning:
## ! NAs introduced by coercion
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 2 remaining warnings.

3 all measures

df %>% 
  ggplot(aes(sample, na_cl_conc_mmolL_TIT))+
  geom_boxplot()+
  geom_point()

df %>% 
  ggplot(aes(sample, na_cl_conc_mmolL_ISE))+
  geom_boxplot()+
  geom_point()
## Warning: Removed 17 rows containing non-finite outside the scale range
## (`stat_boxplot()`).
## Warning: Removed 17 rows containing missing values or values outside the scale range
## (`geom_point()`).

df %>% 
  ggplot(aes(sample, na_cl_conc_mmolL_COND))+
  geom_boxplot()+
  geom_point()
## Warning: Removed 14 rows containing non-finite outside the scale range
## (`stat_boxplot()`).
## Warning: Removed 14 rows containing missing values or values outside the scale range
## (`geom_point()`).

4 only confidence + and ?

df_titr_PQ <- df %>% 
  filter(confidence_titration=="+" |confidence_titration=="?") 

df_titr_PQ%>% 
  ggplot(aes(sample, na_cl_conc_mmolL_TIT))+
  geom_boxplot()+
  geom_point(aes(color = patient))+
  scale_color_manual(values = id_palette)+
  geom_line(aes(group = patient), color="grey")+
  theme_bw()

summary(lmerTest::lmer(na_cl_conc_mmolL_TIT ~ sample + (1|patient), data = df_titr_PQ))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: na_cl_conc_mmolL_TIT ~ sample + (1 | patient)
##    Data: df_titr_PQ
## 
## REML criterion at convergence: 243.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.22336 -0.67439 -0.06745  0.66520  1.74953 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  patient  (Intercept)  328.7   18.13   
##  Residual             1353.9   36.80   
## Number of obs: 32, groups:  patient, 14
## 
## Fixed effects:
##             Estimate Std. Error     df t value Pr(>|t|)    
## (Intercept)    58.00      11.29  18.19   5.136  6.7e-05 ***
## sampleV2      -49.47      25.14  20.33  -1.968   0.0629 .  
## sampleV3      -27.91      29.98  20.34  -0.931   0.3628    
## sampleV4       10.32      30.47  22.98   0.339   0.7380    
## sampleV5       67.15      40.00  15.36   1.679   0.1134    
## sampleV6       42.06      40.00  15.36   1.052   0.3092    
## sampleV7       14.19      28.94  14.86   0.490   0.6310    
## sampleV8      -33.76      24.84  18.43  -1.359   0.1905    
## sampleV9      -12.61      21.85  16.98  -0.577   0.5713    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr) smplV2 smplV3 smplV4 smplV5 smplV6 smplV7 smplV8
## sampleV2 -0.384                                                 
## sampleV3 -0.330  0.243                                          
## sampleV4 -0.334  0.128  0.110                                   
## sampleV5 -0.195  0.078  0.068  0.065                            
## sampleV6 -0.195  0.078  0.068  0.065  0.154                     
## sampleV7 -0.270  0.112  0.099  0.090  0.138  0.138              
## sampleV8 -0.371  0.240  0.238  0.124  0.080  0.080  0.154       
## sampleV9 -0.408  0.197  0.183  0.136  0.140  0.140  0.157  0.194
df_analytics_PQ <- df%>% 
  filter(confidence_anlytics=="+" |confidence_anlytics=="?") 

df_analytics_PQ %>% 
  ggplot(aes(sample, na_cl_conc_mmolL_ISE))+
  geom_boxplot()+
  geom_point(aes(color = patient))+
  scale_color_manual(values = id_palette)+
  geom_line(aes(group = patient), color="grey")+
  theme_classic()

summary(lmerTest::lmer(na_cl_conc_mmolL_ISE ~ sample + (1|patient), data = df_analytics_PQ))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: na_cl_conc_mmolL_ISE ~ sample + (1 | patient)
##    Data: df_analytics_PQ
## 
## REML criterion at convergence: 255.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.5474 -0.3818  0.0000  0.3121  1.5250 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  patient  (Intercept) 1618.8   40.23   
##  Residual              764.1   27.64   
## Number of obs: 33, groups:  patient, 15
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)   68.540     13.056  14.437   5.250 0.000111 ***
## sampleV2     -17.215     21.461  10.911  -0.802 0.439574    
## sampleV3     -16.015     25.430  10.553  -0.630 0.542236    
## sampleV4     -23.599     29.786  20.093  -0.792 0.437445    
## sampleV5      58.592     31.474   8.381   1.862 0.098002 .  
## sampleV6     115.471     31.474   8.381   3.669 0.005833 ** 
## sampleV7      44.852     22.777   8.482   1.969 0.082403 .  
## sampleV8       9.087     20.648  10.123   0.440 0.669122    
## sampleV9       7.768     18.129  10.651   0.428 0.676828    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr) smplV2 smplV3 smplV4 smplV5 smplV6 smplV7 smplV8
## sampleV2 -0.242                                                 
## sampleV3 -0.216  0.383                                          
## sampleV4 -0.245  0.059  0.053                                   
## sampleV5 -0.098  0.051  0.051  0.024                            
## sampleV6 -0.098  0.051  0.051  0.024  0.229                     
## sampleV7 -0.136  0.088  0.090  0.033  0.189  0.189              
## sampleV8 -0.228  0.354  0.376  0.056  0.064  0.064  0.172       
## sampleV9 -0.251  0.200  0.207  0.062  0.166  0.166  0.161  0.204
df_analytics_PQ %>% 
   ggplot(aes(sample, na_cl_conc_mmolL_COND))+
  geom_boxplot()+
  geom_point(aes(color = patient))+
   scale_color_manual(values = id_palette)+
  geom_line(aes(group = patient), color="grey")+
  theme_classic()

summary(lmerTest::lmer(na_cl_conc_mmolL_COND~ sample + (1|patient), data = df_analytics_PQ))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: na_cl_conc_mmolL_COND ~ sample + (1 | patient)
##    Data: df_analytics_PQ
## 
## REML criterion at convergence: 257.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.34176 -0.32617 -0.04848  0.32009  2.05709 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  patient  (Intercept)  824.2   28.71   
##  Residual             1281.3   35.79   
## Number of obs: 33, groups:  patient, 15
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)   87.337     12.317  16.076   7.091 2.49e-06 ***
## sampleV2     -31.971     25.810  15.356  -1.239    0.234    
## sampleV3     -25.899     30.757  14.776  -0.842    0.413    
## sampleV4     -25.262     32.551  23.641  -0.776    0.445    
## sampleV5      56.179     39.866   9.780   1.409    0.190    
## sampleV6      87.273     39.866   9.780   2.189    0.054 .  
## sampleV7      40.877     28.796   9.863   1.420    0.187    
## sampleV8       3.715     25.218  13.320   0.147    0.885    
## sampleV9       4.520     22.054  13.066   0.205    0.841    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr) smplV2 smplV3 smplV4 smplV5 smplV6 smplV7 smplV8
## sampleV2 -0.328                                                 
## sampleV3 -0.286  0.305                                          
## sampleV4 -0.298  0.098  0.085                                   
## sampleV5 -0.147  0.059  0.055  0.044                            
## sampleV6 -0.147  0.059  0.055  0.044  0.194                     
## sampleV7 -0.205  0.091  0.086  0.061  0.163  0.163              
## sampleV8 -0.311  0.286  0.299  0.093  0.066  0.066  0.158       
## sampleV9 -0.341  0.191  0.188  0.102  0.149  0.149  0.152  0.190

5 only confidence +

df_titr_P <- df %>% 
  filter(confidence_titration=="+")

df_titr_P %>% 
  ggplot(aes(sample, na_cl_conc_mmolL_TIT))+
  geom_boxplot()+
  geom_point(aes(color = patient))+
  scale_color_manual(values = id_palette)+
  geom_line(aes(group = patient), color="grey")+
  theme_bw()

summary(lmerTest::lmer(na_cl_conc_mmolL_TIT~ sample + (1|patient), data = df_titr_P)) 
## boundary (singular) fit: see help('isSingular')
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: na_cl_conc_mmolL_TIT ~ sample + (1 | patient)
##    Data: df_titr_P
## 
## REML criterion at convergence: 102.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.3024 -0.6348  0.0000  0.7492  1.3402 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  patient  (Intercept)    0      0.00   
##  Residual             1186     34.43   
## Number of obs: 17, groups:  patient, 11
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)   67.655     11.477  10.000   5.895 0.000152 ***
## sampleV3     -13.749     36.295  10.000  -0.379 0.712750    
## sampleV4      -3.379     26.917  10.000  -0.126 0.902596    
## sampleV5      67.321     36.295  10.000   1.855 0.093302 .  
## sampleV6      42.234     36.295  10.000   1.164 0.271588    
## sampleV7      41.676     36.295  10.000   1.148 0.277590    
## sampleV9      13.356     26.917  10.000   0.496 0.630490    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr) smplV3 smplV4 smplV5 smplV6 smplV7
## sampleV3 -0.316                                   
## sampleV4 -0.426  0.135                            
## sampleV5 -0.316  0.100  0.135                     
## sampleV6 -0.316  0.100  0.135  0.100              
## sampleV7 -0.316  0.100  0.135  0.100  0.100       
## sampleV9 -0.426  0.135  0.182  0.135  0.135  0.135
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
df_ana_P <- df %>% 
  filter(confidence_anlytics=="+") 

df_ana_P%>% 
  ggplot(aes(sample, na_cl_conc_mmolL_ISE))+
  geom_boxplot()+
  geom_point(aes(color = patient))+
  scale_color_manual(values = id_palette)+
  geom_line(aes(group = patient), color="grey")+
  theme_classic()

summary(lmerTest::lmer(na_cl_conc_mmolL_ISE ~ sample + (1|patient), data = df_ana_P)) 
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: na_cl_conc_mmolL_ISE ~ sample + (1 | patient)
##    Data: df_ana_P
## 
## REML criterion at convergence: 244.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.55203 -0.33136 -0.02004  0.15885  1.44578 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  patient  (Intercept) 1728.8   41.58   
##  Residual              717.4   26.79   
## Number of obs: 32, groups:  patient, 15
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  67.9608    13.2157  14.3627   5.142 0.000138 ***
## sampleV2     -0.0378    25.6025  10.3724  -0.001 0.998850    
## sampleV3     -8.0469    25.6025  10.3724  -0.314 0.759522    
## sampleV4    -23.0560    29.4091  17.9571  -0.784 0.443274    
## sampleV5     59.7943    30.5556   7.7321   1.957 0.087323 .  
## sampleV6    116.6731    30.5556   7.7321   3.818 0.005437 ** 
## sampleV7     46.4910    22.1358   7.8446   2.100 0.069592 .  
## sampleV8     15.0906    20.6580   9.8985   0.730 0.482018    
## sampleV9     10.8495    17.7919   9.8658   0.610 0.555776    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr) smplV2 smplV3 smplV4 smplV5 smplV6 smplV7 smplV8
## sampleV2 -0.209                                                 
## sampleV3 -0.209  0.453                                          
## sampleV4 -0.238  0.050  0.050                                   
## sampleV5 -0.094  0.055  0.055  0.022                            
## sampleV6 -0.094  0.055  0.055  0.022  0.232                     
## sampleV7 -0.131  0.099  0.099  0.031  0.191  0.191              
## sampleV8 -0.220  0.418  0.418  0.052  0.068  0.068  0.179       
## sampleV9 -0.243  0.230  0.230  0.058  0.169  0.169  0.166  0.225
df_ana_P %>% 
filter(confidence_anlytics=="+") %>% 
   ggplot(aes(sample, na_cl_conc_mmolL_COND))+
  geom_boxplot()+
  geom_point(aes(color = patient))+
   scale_color_manual(values = id_palette)+
  geom_line(aes(group = patient), color="grey")+
  theme_classic()

summary(lmerTest::lmer(na_cl_conc_mmolL_COND ~ sample + (1|patient), data = df_ana_P)) 
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: na_cl_conc_mmolL_COND ~ sample + (1 | patient)
##    Data: df_ana_P
## 
## REML criterion at convergence: 246.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.33021 -0.35838 -0.01081  0.23264  2.10295 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  patient  (Intercept)  783.2   27.99   
##  Residual             1288.0   35.89   
## Number of obs: 32, groups:  patient, 15
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)   86.965     12.217  14.956   7.119 3.57e-06 ***
## sampleV2     -12.810     31.129  15.292  -0.411   0.6864    
## sampleV3     -20.468     31.129  15.292  -0.658   0.5206    
## sampleV4     -24.991     32.412  22.508  -0.771   0.4487    
## sampleV5      56.466     39.921   8.453   1.414   0.1930    
## sampleV6      87.560     39.921   8.453   2.193   0.0578 .  
## sampleV7      41.367     28.839   8.546   1.434   0.1870    
## sampleV8       7.567     25.465  13.309   0.297   0.7709    
## sampleV9       5.842     22.098  11.801   0.264   0.7960    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr) smplV2 smplV3 smplV4 smplV5 smplV6 smplV7 smplV8
## sampleV2 -0.290                                                 
## sampleV3 -0.290  0.335                                          
## sampleV4 -0.300  0.087  0.087                                   
## sampleV5 -0.150  0.056  0.056  0.045                            
## sampleV6 -0.150  0.056  0.056  0.045  0.192                     
## sampleV7 -0.208  0.089  0.089  0.063  0.162  0.162              
## sampleV8 -0.315  0.311  0.311  0.095  0.068  0.068  0.159       
## sampleV9 -0.345  0.195  0.195  0.104  0.148  0.148  0.153  0.197
df_titr_PQ %>% 
  ggplot(aes(as.numeric(visit), na_cl_conc_mmolL_TIT))+
  geom_point()+
  geom_smooth()
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

df_analytics_PQ %>% 
  ggplot(aes(as.numeric(visit), na_cl_conc_mmolL_ISE))+
  geom_point()+
  geom_smooth()
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

df_analytics_PQ %>% 
  ggplot(aes(as.numeric(visit), na_cl_conc_mmolL_COND))+
  geom_point()+
  geom_smooth()
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

6 only confidence + and ? by patient

df %>% 
  filter(confidence_titration=="+" |confidence_titration=="?") %>% 
  ggplot(aes(sample, na_cl_conc_mmolL_TIT))+
  geom_boxplot()+
  geom_point(aes(color = patient))+
  scale_color_manual(values = id_palette)+
  geom_line(aes(group = patient), color="grey")+
  theme_bw()+
  facet_wrap(~patient)
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?

df %>% 
  filter(confidence_anlytics=="+" |confidence_anlytics=="?") %>% 
  ggplot(aes(sample, na_cl_conc_mmolL_ISE))+
  geom_boxplot()+
  geom_point(aes(color = patient))+
  scale_color_manual(values = id_palette)+
  geom_line(aes(group = patient), color="grey")+
  theme_bw()+
  facet_wrap(~patient)
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?

df %>% 
  filter(confidence_anlytics=="+" |confidence_anlytics=="?") %>% 
  ggplot(aes(sample, na_cl_conc_mmolL_COND))+
  geom_boxplot()+
  geom_point(aes(color = patient))+
  scale_color_manual(values = id_palette)+
  geom_line(aes(group = patient), color="grey")+
  theme_bw()+
  facet_wrap(~patient)
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?

# only confidence + and by patient

df %>% 
  filter(confidence_titration=="+") %>% 
  ggplot(aes(sample, na_cl_conc_mmolL_TIT))+
  geom_boxplot()+
  geom_point(aes(color = patient))+
  scale_color_manual(values = id_palette)+
  geom_line(aes(group = patient), color="grey")+
  theme_bw()+
  facet_wrap(~patient)
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?

df %>% 
  filter(confidence_anlytics=="+") %>% 
  ggplot(aes(sample, na_cl_conc_mmolL_ISE))+
  geom_boxplot()+
  geom_point(aes(color = patient))+
  scale_color_manual(values = id_palette)+
  geom_line(aes(group = patient), color="grey")+
  theme_bw()+
  facet_wrap(~patient)
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?

df %>% 
  filter(confidence_anlytics=="+") %>% 
  ggplot(aes(sample, na_cl_conc_mmolL_COND))+
  geom_boxplot()+
  geom_point(aes(color = patient))+
  scale_color_manual(values = id_palette)+
  geom_line(aes(group = patient), color="grey")+
  theme_bw()+
  facet_wrap(~patient)
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?

7 test correlation of methods

df_fil <- df %>% 
  filter(!is.na(confidence_anlytics)) %>% 
  filter(!is.na(confidence_titration)) %>% 
  filter(confidence_anlytics!="-")

df_fil %>% 
  ggplot(aes(na_cl_conc_mmolL_TIT, na_cl_conc_mmolL_ISE))+
  geom_point()+
  geom_smooth(method = "lm")
## `geom_smooth()` using formula = 'y ~ x'

df_fil %>% 
  ggplot(aes(na_cl_conc_mmolL_TIT, na_cl_conc_mmolL_COND))+
  geom_point()+
  geom_smooth(method = "lm")
## `geom_smooth()` using formula = 'y ~ x'

df_fil %>% 
  ggplot(aes(na_cl_conc_mmolL_ISE, na_cl_conc_mmolL_COND))+
  geom_point()+
  geom_smooth(method = "lm")
## `geom_smooth()` using formula = 'y ~ x'

8 combine NaCl measurements with microbiome and clinical data

ps_sputum_full <- subset_samples(ps_sputum_full, material=="Sputum")
ps_sputum_filtered <- tax_filter(ps_sputum_full, min_prevalence = 0.25) # has 91 taxa
## Proportional min_prevalence given: 0.25 --> min 15/57 samples.
ps_sputum_relab <- transform_sample_counts(ps_sputum_filtered, function(x) x/sum(x))

df_fil<- df_fil %>% 
  mutate(id_visit=paste("IMP", patient, sample, sep = ""))

df_fil_sel <- df_fil %>% 
  select(eppi_n, patient, visit, id_visit, na_cl_conc_mmolL_COND, na_cl_conc_mmolL_ISE, na_cl_conc_mmolL_TIT, confidence_titration, confidence_anlytics)

# extract Staphylococcus abundance from Sputum

sputum_df <- psmelt(ps_sputum_relab)
sputum_df <- sputum_df %>% 
  mutate(Abundance=Abundance*100) # transform relative abundance into percentages

staph_df <- sputum_df%>% 
  filter(Genus=="Staphylococcus")

df_staph <- df_fil_sel %>% 
  left_join(staph_df, by="id_visit")

9 Staphylococcus vs chloride levels in sputum

df_staph %>% 
  ggplot(aes(chlorid_mmol_cl_l, na_cl_conc_mmolL_COND))+
  geom_point()+
  geom_smooth()
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## Warning: Removed 10 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_point()`).

  library(ggtext)
## Warning: package 'ggtext' was built under R version 4.1.2
visit_sum_palette <- c("black", "#C3BC3FFF","#6388B4FF","#BB7693FF","#55AD89FF") 
x2labels <- c("0", "3", "6-12", "15-18", "21-24")

df_staph %>% 
ggplot(aes(y = Abundance+ 0.01, x=chlorid_mmol_cl_l)) +
  geom_point(aes(color = visit_sum), size = 4) +
  geom_smooth(method = "lm", color = "black", se = T) +
  scale_color_manual(values = visit_sum_palette, labels = x2labels) +
  scale_y_log10(
    breaks = c(0.01, 0.1, 1, 10, 100),
    labels = c("0", "0.1", "1", "10", "100"), limits = c(0.005, 100)) +
  xlab("Sweat chloride mmol/l") +
  ylab("*Staphylococcus* (%) in Sputum") +
  theme_classic() +
  theme(
    axis.title.y  = ggtext::element_markdown(size = 12),
    axis.title.x  = ggtext::element_markdown(size = 12),
    legend.title = element_text(size = 12),
    text = element_text(size = 12),
    legend.position = "bottom",
    axis.text = element_text(size = 12)
  ) +
  annotate("text", x = 90, y = 200, label = "fdr=0.015, Ds=0.57", size = 5) +
  geom_vline(xintercept = 30, linetype = 'dotted', col = 'red') +
  geom_vline(xintercept = 60, linetype = 'dotted', col = 'red') +
  guides(color = guide_legend(title = "Months from\ntreatment start", title.position = "left", ncol = 3)) 
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 10 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 10 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 5 rows containing missing values or values outside the scale range
## (`geom_smooth()`).
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_text()`).

df_staph %>% 
ggplot(aes(y = Abundance+ 0.01, x=na_cl_conc_mmolL_COND)) +
  geom_point(aes(color = visit_sum), size = 4) +
  geom_smooth(method = "lm", color = "black", se = T) +
  scale_color_manual(values = visit_sum_palette, labels = x2labels) +
  scale_y_log10(
    breaks = c(0.01, 0.1, 1, 10, 100),
    labels = c("0", "0.1", "1", "10", "100"), limits = c(0.005, 100)) +
  xlab("Chloride mmol/l (conductivity) in Sputum") +
  ylab("*Staphylococcus* (%) in Sputum") +
  theme_classic() +
  theme(
    axis.title.y  = ggtext::element_markdown(size = 12),
    axis.title.x  = ggtext::element_markdown(size = 12),
    legend.title = element_text(size = 12),
    text = element_text(size = 12),
    legend.position = "bottom",
    axis.text = element_text(size = 12)
  ) +
  #annotate("text", x = 90, y = 200, label = "fdr=0.015, Ds=0.57", size = 5) +
  #geom_vline(xintercept = 30, linetype = 'dotted', col = 'red') +
  #geom_vline(xintercept = 60, linetype = 'dotted', col = 'red') +
  guides(color = guide_legend(title = "Months from\ntreatment start", title.position = "left", ncol = 3)) 
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 3 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).

df_staph %>% 
ggplot(aes(y = Abundance+ 0.01, x=na_cl_conc_mmolL_ISE)) +
  geom_point(aes(color = visit_sum), size = 4) +
  geom_smooth(method = "lm", color = "black", se = T) +
  scale_color_manual(values = visit_sum_palette, labels = x2labels) +
  scale_y_log10(
    breaks = c(0.01, 0.1, 1, 10, 100),
    labels = c("0", "0.1", "1", "10", "100"), limits = c(0.005, 100)) +
  xlab("Chloride mmol/l (ISE) in Sputum") +
  ylab("*Staphylococcus* (%) in Sputum") +
  theme_classic() +
  theme(
    axis.title.y  = ggtext::element_markdown(size = 12),
    axis.title.x  = ggtext::element_markdown(size = 12),
    legend.title = element_text(size = 12),
    text = element_text(size = 12),
    legend.position = "bottom",
    axis.text = element_text(size = 12)
  ) +
  #annotate("text", x = 90, y = 200, label = "fdr=0.015, Ds=0.57", size = 5) +
  #geom_vline(xintercept = 30, linetype = 'dotted', col = 'red') +
  #geom_vline(xintercept = 60, linetype = 'dotted', col = 'red') +
  guides(color = guide_legend(title = "Months from\ntreatment start", title.position = "left", ncol = 3)) 
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 3 rows containing non-finite outside the scale range (`stat_smooth()`).
## Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).

df_staph %>% 
ggplot(aes(y = Abundance+ 0.01, x=na_cl_conc_mmolL_TIT)) +
  geom_point(aes(color = visit_sum), size = 4) +
  geom_smooth(method = "lm", color = "black", se = T) +
  scale_color_manual(values = visit_sum_palette, labels = x2labels) +
  scale_y_log10(
    breaks = c(0.01, 0.1, 1, 10, 100),
    labels = c("0", "0.1", "1", "10", "100"), limits = c(0.005, 100)) +
  xlab("Chloride mmol/l (titration) in Sputum") +
  ylab("*Staphylococcus* (%) in Sputum") +
  theme_classic() +
  theme(
    axis.title.y  = ggtext::element_markdown(size = 12),
    axis.title.x  = ggtext::element_markdown(size = 12),
    legend.title = element_text(size = 12),
    text = element_text(size = 12),
    legend.position = "bottom",
    axis.text = element_text(size = 12)
  ) +
  #annotate("text", x = 90, y = 200, label = "fdr=0.015, Ds=0.57", size = 5) +
  #geom_vline(xintercept = 30, linetype = 'dotted', col = 'red') +
  #geom_vline(xintercept = 60, linetype = 'dotted', col = 'red') +
  guides(color = guide_legend(title = "Months from\ntreatment start", title.position = "left", ncol = 3)) 
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 3 rows containing non-finite outside the scale range (`stat_smooth()`).
## Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).